```r
require(phyloseq)
require(tidyverse)
require(phyloseq)
require(reshape2)
require(dplyr)
require(ggplot2)
require(microbiome)
<!-- rnb-source-end -->
<!-- rnb-output-begin eyJkYXRhIjoiTG9hZGluZyByZXF1aXJlZCBwYWNrYWdlOiBtaWNyb2Jpb21lXG5XYXJuaW5nIGluIGxpYnJhcnkocGFja2FnZSwgbGliLmxvYyA9IGxpYi5sb2MsIGNoYXJhY3Rlci5vbmx5ID0gVFJVRSwgbG9naWNhbC5yZXR1cm4gPSBUUlVFLCAgOlxuICB0aGVyZSBpcyBubyBwYWNrYWdlIGNhbGxlZCDigJhtaWNyb2Jpb21l4oCZXG4ifQ== -->
Loading required package: microbiome Warning in library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE, : there is no package called ‘microbiome’
<!-- rnb-output-end -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucmVxdWlyZSh2ZWdhbilcbmBgYFxuYGBgIn0= -->
```r
```r
require(vegan)
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Load data
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxucHNfZG1uIDwtIHJlYWRSRFMoXCIvVXNlcnMvZ29yZG9uY3VzdGVyL0Rlc2t0b3AvR2l0X1Byb2plY3RzL0hlcmJpY2lkZV9NaWNyb2Jlc19QVDEvZGF0YS9QaHlsb3NlcU9iamVjdHMvSVRTL0RNTl9lc3RzX0lUUy5SZGF0YVwiKVxuYGBgIn0= -->
```r
ps_dmn <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/ITS/DMN_ests_ITS.Rdata")
Warning message:
In rm(list = objects, envir = env) : object '*tmp*' not found
sample_data(ps_dmn)$Herbicide <- factor(sample_data(ps_dmn)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_dmn)$herb_time<-paste(sample_data(ps_dmn)$Herbicide, sample_data(ps_dmn)$Time, sep = "_")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_dmn)$Mode<-sample_data(ps_dmn)$Herbicide
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")
sample_data(ps_dmn)$Mode<- as.factor(values[match(sample_data(ps_dmn)$Mode, index)])
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
sample_data(ps_dmn)$Herbicide <- as.factor(values[match(sample_data(ps_dmn)$Herbicide, index)])
ps_rare <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/ITS/HerbPt1_rare_ITS.Rdata")
sample_data(ps_rare)$Herbicide <- factor(sample_data(ps_rare)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_rare)$herb_time<-paste(sample_data(ps_rare)$Herbicide, sample_data(ps_rare)$Time, sep = "_")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")
sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
sample_data(ps_rare)$Herbicide <- as.factor(values[match(sample_data(ps_rare)$Herbicide, index)])
ps_trans <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/ITS/HerbPt1_hel_trans_ITS.Rdata")
sample_data(ps_trans)$Herbicide <- factor(sample_data(ps_trans)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_trans)$herb_time<-paste(sample_data(ps_trans)$Herbicide, sample_data(ps_trans)$Time, sep = "_")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_trans)$Mode<-sample_data(ps_trans)$Herbicide
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")
sample_data(ps_trans)$Mode<- as.factor(values[match(sample_data(ps_trans)$Mode, index)])
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
sample_data(ps_trans)$Herbicide <- as.factor(values[match(sample_data(ps_trans)$Herbicide, index)])
Remove samples that are outlines or under sequenced.
ps_dmn <- subset_samples(ps_dmn, sample_names(ps_dmn) != "G009SG")
ps_dmn <- subset_samples(ps_dmn, sample_names(ps_dmn) != "G095SG")
ps_dmn <- subset_samples(ps_dmn, sample_names(ps_dmn) != "G123SG")
ps_dmn <- subset_samples(ps_dmn, sample_names(ps_dmn) != "G129SG")
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G009SG")
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G095SG")
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G123SG")
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G129SG")
ps_rare_sub<-prune_taxa(taxa_sums(ps_rare) > 2, ps_rare)
ordinations and adonis testing with three separate objects (i.e., dmn, rarefied, transformed). Rare taxa are removed from rarefied and transformed to successfully ordinate. At this point, the transformed data will not ordinate. This section is full dataset ordinations.
ord_dmn<-ordinate(physeq = ps_dmn, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit = 1000)
Run 0 stress 0.1637061
Run 1 stress 0.1638069
... Procrustes: rmse 0.04511874 max resid 0.2959219
Run 2 stress 0.1639952
... Procrustes: rmse 0.04572109 max resid 0.2996605
Run 3 stress 0.1636993
... New best solution
... Procrustes: rmse 0.003124 max resid 0.03552253
Run 4 stress 0.1636321
... New best solution
... Procrustes: rmse 0.01768586 max resid 0.1058303
Run 5 stress 0.1634893
... New best solution
... Procrustes: rmse 0.01749345 max resid 0.1018927
Run 6 stress 0.1636166
... Procrustes: rmse 0.04806717 max resid 0.3003009
Run 7 stress 0.1633705
... New best solution
... Procrustes: rmse 0.04470724 max resid 0.2996104
Run 8 stress 0.1631225
... New best solution
... Procrustes: rmse 0.01730151 max resid 0.1523431
Run 9 stress 0.1636991
Run 10 stress 0.1640706
Run 11 stress 0.1635635
... Procrustes: rmse 0.01665523 max resid 0.203607
Run 12 stress 0.163487
... Procrustes: rmse 0.04938977 max resid 0.2995861
Run 13 stress 0.1642535
Run 14 stress 0.1633718
... Procrustes: rmse 0.01725707 max resid 0.1522712
Run 15 stress 0.1640259
Run 16 stress 0.1646096
Run 17 stress 0.1634884
... Procrustes: rmse 0.04947432 max resid 0.2997511
Run 18 stress 0.1631224
... New best solution
... Procrustes: rmse 0.001164614 max resid 0.007901557
... Similar to previous best
Run 19 stress 0.1631223
... New best solution
... Procrustes: rmse 0.000894574 max resid 0.007191738
... Similar to previous best
Run 20 stress 0.1632532
... Procrustes: rmse 0.007653146 max resid 0.08811271
*** Solution reached
ord_rare<-ordinate(physeq = ps_rare_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit = 1000)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1837394
Run 1 stress 0.1837536
... Procrustes: rmse 0.004222286 max resid 0.03854408
Run 2 stress 0.1859057
Run 3 stress 0.1837379
... New best solution
... Procrustes: rmse 0.0006553106 max resid 0.004505322
... Similar to previous best
Run 4 stress 0.1856626
Run 5 stress 0.1861936
Run 6 stress 0.1862078
Run 7 stress 0.1859133
Run 8 stress 0.1837401
... Procrustes: rmse 0.0007329372 max resid 0.004394223
... Similar to previous best
Run 9 stress 0.1837556
... Procrustes: rmse 0.004232884 max resid 0.03808415
Run 10 stress 0.1837291
... New best solution
... Procrustes: rmse 0.004668855 max resid 0.03839136
Run 11 stress 0.1862108
Run 12 stress 0.1856156
Run 13 stress 0.1837321
... Procrustes: rmse 0.001869121 max resid 0.0172899
Run 14 stress 0.185933
Run 15 stress 0.1859134
Run 16 stress 0.1837366
... Procrustes: rmse 0.003934702 max resid 0.03800597
Run 17 stress 0.183736
... Procrustes: rmse 0.004072096 max resid 0.03782902
Run 18 stress 0.1862073
Run 19 stress 0.1837283
... New best solution
... Procrustes: rmse 0.003314344 max resid 0.02283046
Run 20 stress 0.1861278
Run 21 stress 0.1837368
... Procrustes: rmse 0.004202959 max resid 0.03676119
Run 22 stress 0.1862196
Run 23 stress 0.1837526
... Procrustes: rmse 0.002014267 max resid 0.0209574
Run 24 stress 0.1837541
... Procrustes: rmse 0.001993472 max resid 0.02133405
Run 25 stress 0.183726
... New best solution
... Procrustes: rmse 0.0003499443 max resid 0.002495313
... Similar to previous best
*** Solution reached
ps_trans_sub<-prune_taxa(taxa_sums(ps_trans) > 0.01, ps_trans)
ord_transformed<-ordinate(physeq = ps_trans_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit = 1000)
Run 0 stress 0.1050141
Run 1 stress 0.103327
... New best solution
... Procrustes: rmse 0.03564627 max resid 0.1241836
Run 2 stress 0.1033284
... Procrustes: rmse 0.0004026237 max resid 0.003428707
... Similar to previous best
Run 3 stress 0.1033322
... Procrustes: rmse 0.000709747 max resid 0.004562807
... Similar to previous best
Run 4 stress 0.1034339
... Procrustes: rmse 0.00953169 max resid 0.05735493
Run 5 stress 0.103324
... New best solution
... Procrustes: rmse 0.002440901 max resid 0.0170954
Run 6 stress 0.1034354
... Procrustes: rmse 0.007967689 max resid 0.05434794
Run 7 stress 0.1033282
... Procrustes: rmse 0.00265007 max resid 0.01838042
Run 8 stress 0.1034361
... Procrustes: rmse 0.007981759 max resid 0.05521913
Run 9 stress 0.1034359
... Procrustes: rmse 0.008040017 max resid 0.05456366
Run 10 stress 0.1033249
... Procrustes: rmse 0.002039597 max resid 0.01573552
Run 11 stress 0.1050151
Run 12 stress 0.1034358
... Procrustes: rmse 0.00800128 max resid 0.05441324
Run 13 stress 0.1050135
Run 14 stress 0.1068755
Run 15 stress 0.1050132
Run 16 stress 0.1033314
... Procrustes: rmse 0.003020749 max resid 0.02047583
Run 17 stress 0.1050126
Run 18 stress 0.103327
... Procrustes: rmse 0.002759821 max resid 0.02293387
Run 19 stress 0.1033299
... Procrustes: rmse 0.002844598 max resid 0.01951695
Run 20 stress 0.103436
... Procrustes: rmse 0.008047048 max resid 0.05456818
Run 21 stress 0.1050133
Run 22 stress 0.1050145
Run 23 stress 0.1034335
... Procrustes: rmse 0.007696545 max resid 0.05320767
Run 24 stress 0.1050135
Run 25 stress 0.1034323
... Procrustes: rmse 0.007414607 max resid 0.05157759
Run 26 stress 0.1033298
... Procrustes: rmse 0.002836647 max resid 0.01947178
Run 27 stress 0.1050129
Run 28 stress 0.1050133
Run 29 stress 0.1050146
Run 30 stress 0.1033263
... Procrustes: rmse 0.002509226 max resid 0.02018023
Run 31 stress 0.103322
... New best solution
... Procrustes: rmse 0.00036451 max resid 0.002650074
... Similar to previous best
*** Solution reached
create alphadiversity tables
alpha_div_md$Herbicide
[1] Dicamba Dicamba Dicamba Handweeded Handweeded
[6] Handweeded Non-Treated Non-Treated Atrazine-Mesotrione Atrazine-Mesotrione
[11] Glyphosate Glyphosate Glyphosate Glyphosate Glyphosate
[16] Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded Handweeded
[21] Dicamba Dicamba Dicamba Non-Treated Non-Treated
[26] Non-Treated Non-Treated Non-Treated Non-Treated Glyphosate
[31] Glyphosate Glyphosate Dicamba Dicamba Dicamba
[36] Atrazine-Mesotrione Atrazine-Mesotrione Handweeded Handweeded Handweeded
[41] Glyphosate Glyphosate Glyphosate Handweeded Handweeded
[46] Non-Treated Non-Treated Non-Treated Dicamba Dicamba
[51] Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Dicamba Dicamba
[56] Handweeded Handweeded Handweeded Non-Treated Non-Treated
[61] Non-Treated Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Glyphosate
[66] Glyphosate Glyphosate Glyphosate Glyphosate Glyphosate
[71] Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded Handweeded
[76] Handweeded Dicamba Dicamba Dicamba Non-Treated
[81] Non-Treated Non-Treated Non-Treated Glyphosate Glyphosate
[86] Dicamba Dicamba Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione
[91] Handweeded Handweeded Handweeded Glyphosate Glyphosate
[96] Glyphosate Handweeded Handweeded Non-Treated Non-Treated
[101] Non-Treated Dicamba Dicamba Atrazine-Mesotrione Dicamba
[106] Dicamba Handweeded Handweeded Handweeded Non-Treated
[111] Non-Treated Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Glyphosate
[116] Glyphosate Glyphosate Glyphosate Glyphosate Atrazine-Mesotrione
[121] Atrazine-Mesotrione Atrazine-Mesotrione Handweeded Handweeded Handweeded
[126] Dicamba Dicamba Dicamba Non-Treated Non-Treated
[131] Non-Treated Non-Treated Non-Treated Glyphosate Glyphosate
[136] Glyphosate Dicamba Dicamba Dicamba Atrazine-Mesotrione
[141] Atrazine-Mesotrione Atrazine-Mesotrione Handweeded Handweeded Glyphosate
[146] Glyphosate Glyphosate Handweeded Handweeded Handweeded
[151] Non-Treated Non-Treated Non-Treated Dicamba Atrazine-Mesotrione
[156] Atrazine-Mesotrione Atrazine-Mesotrione
Levels: Non-Treated Handweeded Atrazine-Mesotrione Dicamba Glyphosate
Shannon Div plots - no significant differences among herbicide treatments at any of the three time points
ggplot(data = alpha_div_md, aes(Herbicide, Shannon, color= Herbicide)) + facet_grid(. ~ Time) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1) )
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_Shannon.pdf")
Saving 7.29 x 4.51 in image
aov_t1<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T1",])
plot(aov_t1$residuals)
summary(aov_t1)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 0.389 0.09720 1.26 0.299
Residuals 48 3.704 0.07716
aov_t2<-aov(Shannon~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T2",])
plot(aov_t2$residuals)
summary(aov_t2)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 0.1234 0.03086 0.508 0.73
Residuals 46 2.7936 0.06073
aov_t3<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3",])
plot(aov_t3$residuals)
summary(aov_t3)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 0.5376 0.13440 2.551 0.051 .
Residuals 48 2.5287 0.05268
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(aov_t3, "Herbicide")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3", ])
$Herbicide
diff lwr upr p adj
Handweeded-Non-Treated 0.086844630 -0.19738244 0.3710717 0.9078666
Atrazine-Mesotrione-Non-Treated 0.129708505 -0.14882205 0.4082391 0.6803851
Dicamba-Non-Treated -0.002861275 -0.30174872 0.2960262 0.9999999
Glyphosate-Non-Treated 0.273363070 -0.01086400 0.5575901 0.0647237
Atrazine-Mesotrione-Handweeded 0.042863876 -0.22867317 0.3144009 0.9914461
Dicamba-Handweeded -0.089705904 -0.38208716 0.2026754 0.9066040
Glyphosate-Handweeded 0.186518440 -0.09085878 0.4638957 0.3282566
Dicamba-Atrazine-Mesotrione -0.132569780 -0.41941651 0.1542769 0.6864852
Glyphosate-Atrazine-Mesotrione 0.143654565 -0.12788248 0.4151916 0.5679081
Glyphosate-Dicamba 0.276224345 -0.01615691 0.5686056 0.0724032
Adonis testing of herbicide treatments by time point
ps_adonis<-function(physeq){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
print(anova(vegan::betadisper(physeq_dist, md_tab$Herbicide)))
vegan::adonis(physeq_dist ~ Herbicide * Time + Total_Weed_Veg , data = md_tab, permutations = 1000)
}
#ps_adonis(ps_rare_sub)
#remove one sample with no vegetation measurement.
ps_rare_sub_57<-subset_samples(ps_rare_sub, sample_names(ps_rare_sub) != "G065SG")
ps_adonis(ps_rare_sub_57)
#ps_adonis(ps_trans_sub)
ps_dmn_57<-subset_samples(ps_dmn, sample_names(ps_dmn) != "G065SG")
#ps_adonis(ps_dmn)
ps_adonis(ps_dmn_57)
Ordination plots DMN by time point
ord_t1_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Run 0 stress 0.1335559
Run 1 stress 0.1309245
... New best solution
... Procrustes: rmse 0.06976384 max resid 0.2060325
Run 2 stress 0.1315098
Run 3 stress 0.132755
Run 4 stress 0.134753
Run 5 stress 0.1308758
... New best solution
... Procrustes: rmse 0.0538814 max resid 0.2344595
Run 6 stress 0.1331698
Run 7 stress 0.1323904
Run 8 stress 0.1346675
Run 9 stress 0.1350022
Run 10 stress 0.1322358
Run 11 stress 0.1334733
Run 12 stress 0.1315248
Run 13 stress 0.1331407
Run 14 stress 0.1327605
Run 15 stress 0.1332143
Run 16 stress 0.1369199
Run 17 stress 0.1339829
Run 18 stress 0.1331858
Run 19 stress 0.1360529
Run 20 stress 0.1352334
Run 21 stress 0.1308759
... Procrustes: rmse 0.001713885 max resid 0.006848176
... Similar to previous best
*** Solution reached
T1_dmn<-ggordiplots::gg_ordiplot(ord = ord_t1_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_dmn_T1.pdf")
Saving 7.29 x 4.51 in image
ord_t2_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Run 0 stress 0.1325262
Run 1 stress 0.1325269
... Procrustes: rmse 0.0001074309 max resid 0.000408716
... Similar to previous best
Run 2 stress 0.1325262
... New best solution
... Procrustes: rmse 0.001015874 max resid 0.005488789
... Similar to previous best
Run 3 stress 0.1327962
... Procrustes: rmse 0.01038273 max resid 0.06195753
Run 4 stress 0.1338926
Run 5 stress 0.1392078
Run 6 stress 0.1338929
Run 7 stress 0.1325264
... Procrustes: rmse 0.001066291 max resid 0.005401677
... Similar to previous best
Run 8 stress 0.1325261
... New best solution
... Procrustes: rmse 0.0001868934 max resid 0.001066974
... Similar to previous best
Run 9 stress 0.1334238
Run 10 stress 0.1325259
... New best solution
... Procrustes: rmse 0.0007505591 max resid 0.004597761
... Similar to previous best
Run 11 stress 0.132925
... Procrustes: rmse 0.01264762 max resid 0.07728014
Run 12 stress 0.1438152
Run 13 stress 0.132526
... Procrustes: rmse 0.0003451055 max resid 0.001743723
... Similar to previous best
Run 14 stress 0.1325259
... Procrustes: rmse 0.0003370863 max resid 0.00167214
... Similar to previous best
Run 15 stress 0.1334866
Run 16 stress 0.1325266
... Procrustes: rmse 0.0009514394 max resid 0.005508805
... Similar to previous best
Run 17 stress 0.1325265
... Procrustes: rmse 0.0008955128 max resid 0.005161325
... Similar to previous best
Run 18 stress 0.1449845
Run 19 stress 0.1342762
Run 20 stress 0.1325263
... Procrustes: rmse 0.0003284059 max resid 0.001374251
... Similar to previous best
*** Solution reached
T2_dmn<-ggordiplots::gg_ordiplot(ord = ord_t2_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_dmn_T2.pdf")
Saving 7.29 x 4.51 in image
#this time point needs to be checked out. The ordination is not working properly.
ord_t3_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Run 0 stress 0.1451733
Run 1 stress 0.1459287
Run 2 stress 0.1510555
Run 3 stress 0.1501906
Run 4 stress 0.1503806
Run 5 stress 0.1459215
Run 6 stress 0.1528692
Run 7 stress 0.1451733
... Procrustes: rmse 0.00281782 max resid 0.01088482
Run 8 stress 0.1502075
Run 9 stress 0.1507344
Run 10 stress 0.1478567
Run 11 stress 0.1490906
Run 12 stress 0.1451724
... New best solution
... Procrustes: rmse 0.0002730257 max resid 0.001361369
... Similar to previous best
Run 13 stress 0.1507808
Run 14 stress 0.146477
Run 15 stress 0.1510089
Run 16 stress 0.1507665
Run 17 stress 0.1518463
Run 18 stress 0.1459425
Run 19 stress 0.1451714
... New best solution
... Procrustes: rmse 0.001928599 max resid 0.006972218
... Similar to previous best
Run 20 stress 0.1459215
*** Solution reached
T3_dmn<-ggordiplots::gg_ordiplot(ord = ord_t3_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_dmn_T3.pdf")
Saving 7.29 x 4.51 in image
Ordination plots on rarefied data by time point.
ord_t1_rare<-ordinate(physeq = subset_samples(ps_rare_sub, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.15407
Run 1 stress 0.1547275
Run 2 stress 0.1560278
Run 3 stress 0.154304
... Procrustes: rmse 0.01759432 max resid 0.08547907
Run 4 stress 0.1552432
Run 5 stress 0.1559291
Run 6 stress 0.1560286
Run 7 stress 0.1560271
Run 8 stress 0.1540703
... Procrustes: rmse 0.0005180512 max resid 0.002163285
... Similar to previous best
Run 9 stress 0.1548105
Run 10 stress 0.1540705
... Procrustes: rmse 0.0005653613 max resid 0.002608133
... Similar to previous best
Run 11 stress 0.1581519
Run 12 stress 0.1559287
Run 13 stress 0.1559305
Run 14 stress 0.1595564
Run 15 stress 0.1540698
... New best solution
... Procrustes: rmse 0.0002352417 max resid 0.000845193
... Similar to previous best
Run 16 stress 0.1585595
Run 17 stress 0.1540701
... Procrustes: rmse 0.0004692716 max resid 0.002323505
... Similar to previous best
Run 18 stress 0.1541024
... Procrustes: rmse 0.004124925 max resid 0.02278237
Run 19 stress 0.1543293
... Procrustes: rmse 0.008122859 max resid 0.02901827
Run 20 stress 0.158152
*** Solution reached
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare$plot + theme_classic()
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare_plot<-T1_rare$plot + theme_classic() + xlim(-0.5, 0.5) + ylim(-0.5, 0.5) + guides(color=guide_legend("Treatment")) + xlab("NMDS 1") + ylab("NMDS 2")
T1_rare_plot
Warning: Removed 3 rows containing missing values (geom_point).
library(cowplot)
my_legend <- get_legend(T1_rare_plot)
Warning: Removed 3 rows containing missing values (geom_point).
library(ggpubr)
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordinationlegend.pdf")
Saving 7.29 x 4.51 in image
T1_rare_plot<-T1_rare_plot + theme(legend.position = "none")
T1_rare_plot
Warning: Removed 3 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T1.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 3 rows containing missing values (geom_point).
ord_t2_rare<-ordinate(physeq = subset_samples(ps_rare_sub, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1716481
Run 1 stress 0.1712224
... New best solution
... Procrustes: rmse 0.07764546 max resid 0.2246421
Run 2 stress 0.1706262
... New best solution
... Procrustes: rmse 0.07905606 max resid 0.2661682
Run 3 stress 0.1713448
Run 4 stress 0.1709004
... Procrustes: rmse 0.05656489 max resid 0.2944112
Run 5 stress 0.1702319
... New best solution
... Procrustes: rmse 0.06383277 max resid 0.2810499
Run 6 stress 0.1706227
... Procrustes: rmse 0.06377346 max resid 0.2806419
Run 7 stress 0.1715869
Run 8 stress 0.1713431
Run 9 stress 0.1760045
Run 10 stress 0.1710628
Run 11 stress 0.1713372
Run 12 stress 0.1709006
Run 13 stress 0.1730992
Run 14 stress 0.1721699
Run 15 stress 0.1720581
Run 16 stress 0.1702893
... Procrustes: rmse 0.05501828 max resid 0.2833975
Run 17 stress 0.1724195
Run 18 stress 0.1707193
... Procrustes: rmse 0.01192387 max resid 0.06232137
Run 19 stress 0.170143
... New best solution
... Procrustes: rmse 0.006859545 max resid 0.02883077
Run 20 stress 0.172169
Run 21 stress 0.1703025
... Procrustes: rmse 0.05780913 max resid 0.2825995
Run 22 stress 0.1709835
Run 23 stress 0.1702878
... Procrustes: rmse 0.05725869 max resid 0.2827813
Run 24 stress 0.1706262
... Procrustes: rmse 0.06561008 max resid 0.2790836
Run 25 stress 0.1710053
Run 26 stress 0.1715864
Run 27 stress 0.1770233
Run 28 stress 0.1728199
Run 29 stress 0.1702891
... Procrustes: rmse 0.0574152 max resid 0.2826642
Run 30 stress 0.1736522
Run 31 stress 0.1733041
Run 32 stress 0.1712221
Run 33 stress 0.1720533
Run 34 stress 0.1733034
Run 35 stress 0.1727863
Run 36 stress 0.1762042
Run 37 stress 0.1767485
Run 38 stress 0.1710044
Run 39 stress 0.1703045
... Procrustes: rmse 0.05793853 max resid 0.2825322
Run 40 stress 0.1701603
... Procrustes: rmse 0.00243186 max resid 0.007204675
... Similar to previous best
*** Solution reached
T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_rare_plot<-T2_rare$plot + theme_classic() + xlim(-0.5, 0.5) + ylim(-0.5, 0.5) + theme(legend.position = "none") + xlab("NMDS 1") + ylab("NMDS 2")
T2_rare_plot
Warning: Removed 4 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T2.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 4 rows containing missing values (geom_point).
ord_t3_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1685562
Run 1 stress 0.1748996
Run 2 stress 0.1726844
Run 3 stress 0.1728187
Run 4 stress 0.1708916
Run 5 stress 0.1723652
Run 6 stress 0.1708747
Run 7 stress 0.1718667
Run 8 stress 0.1686699
... Procrustes: rmse 0.009687857 max resid 0.05462871
Run 9 stress 0.1684508
... New best solution
... Procrustes: rmse 0.0119619 max resid 0.06653228
Run 10 stress 0.1772936
Run 11 stress 0.1687617
... Procrustes: rmse 0.02053262 max resid 0.105457
Run 12 stress 0.1686702
... Procrustes: rmse 0.0151362 max resid 0.06254477
Run 13 stress 0.1686729
... Procrustes: rmse 0.02086917 max resid 0.1253609
Run 14 stress 0.1792488
Run 15 stress 0.1694632
Run 16 stress 0.1696437
Run 17 stress 0.170417
Run 18 stress 0.1823142
Run 19 stress 0.1687895
... Procrustes: rmse 0.02300242 max resid 0.1239071
Run 20 stress 0.1687873
... Procrustes: rmse 0.02263201 max resid 0.1209943
Run 21 stress 0.1723288
Run 22 stress 0.1697892
Run 23 stress 0.1694646
Run 24 stress 0.1686697
... Procrustes: rmse 0.01566615 max resid 0.06646188
Run 25 stress 0.1720066
Run 26 stress 0.1708448
Run 27 stress 0.172342
Run 28 stress 0.1708482
Run 29 stress 0.1686102
... Procrustes: rmse 0.01107215 max resid 0.06307271
Run 30 stress 0.1686539
... Procrustes: rmse 0.01829401 max resid 0.1088791
Run 31 stress 0.1725083
Run 32 stress 0.1729314
Run 33 stress 0.1823118
Run 34 stress 0.171603
Run 35 stress 0.1684516
... Procrustes: rmse 0.000333685 max resid 0.001749564
... Similar to previous best
*** Solution reached
T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_rare_plot<-T3_rare$plot + theme_classic() + xlim(-0.5, 0.5) + ylim(-0.5, 0.5) + theme(legend.position = "none") + xlab("NMDS 1") + ylab("NMDS 2")
T3_rare_plot
Warning: Removed 4 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T3.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 4 rows containing missing values (geom_point).
library(ggpubr)
ggarrange(T1_rare_plot, T2_rare_plot, T3_rare_plot, ncol = 1)
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 4 rows containing missing values (geom_point).
Warning: Removed 4 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_ordination.pdf", width = 5, height = 10)
CAP ordination plots rarefied
t1_dist <- distance(subset_samples(ps_rare, Time=="T1"), method="bray") #get wUnifrac and save
t1_table<-as.matrix(dist(t1_dist)) #transform wUnifrac index
ord_t1_rare_cap <- capscale(t1_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T1"))))
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T1_cap.pdf")
t2_dist <- distance(subset_samples(ps_rare, Time=="T2"), method="bray") #get wUnifrac and save
t2_table<-as.matrix(dist(t2_dist)) #transform wUnifrac index
ord_t2_rare_cap <- capscale(t2_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T2"))))
T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T2_cap.pdf")
#G166SG identified as outlier based on plots with it included. Removed to create plot.
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
t3_dist <- distance(subset_samples(ps_rare, Time=="T3"), method="bray") #get wUnifrac and save
t3_table<-as.matrix(dist(t3_dist)) #transform wUnifrac index
ord_t3_rare_cap <- capscale(t3_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T3"))))
T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T3_cap.pdf")
Pairwise adonis testing no needed becasue of insignificant gloabal test.
ps_pairwiseadonis<-function(physeq){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
pairwiseAdonis::pairwise.adonis(x = physeq_dist, factors = md_tab$Herbicide, p.adjust.m = "none", perm = 1000)
}
ps_t1<-subset_samples(ps_rare_sub, Time == "T1")
ps_t1<-prune_taxa(taxa_sums(ps_t1) > 2, ps_t1)
ps_t2<-subset_samples(ps_rare_sub, Time == "T2")
ps_t2<-prune_taxa(taxa_sums(ps_t2) > 2, ps_t2)
ps_t3<-subset_samples(ps_rare_sub, Time == "T3")
ps_t3<-prune_taxa(taxa_sums(ps_t3) > 2, ps_t3)
ps_pairwiseadonis(ps_t1)
ps_pairwiseadonis(ps_t2)
ps_pairwiseadonis(ps_t3)
Pairwise betadispr by treatment, time and mode
ps_betadispr<-function(physeq, groupingvar = "Groupingvar"){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
mod<-vegan::betadisper(physeq_dist, md_tab[,groupingvar])
## Perform test
print(anova(mod))
## Permutation test for F
pmod <- vegan::permutest(mod, permutations = 1000, pairwise = TRUE)
print(pmod)
print(boxplot(mod))
}
permute test of dispersion
ps_betadispr(subset_samples(ps_rare_sub, Time == "T1"), groupingvar = "Mode")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T2"), groupingvar = "Mode")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T3"), groupingvar = "Mode")
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Chemical"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Non-Treated"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Hand"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Glyphosate"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Atrazine-Mesotrione"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Dicamba"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Handweeded"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Non-Treated"), groupingvar = "Time")
ps_betadispr(ps_rare_sub, groupingvar = "Time")
Box and whisker plots of distance within group distances
#remotes::install_github("antonioggsousa/micrUBIfuns")
library(micrUBIfuns)
T1_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Herbicide")
T1_beta_plot <- T1_beta$plot
T1_beta_plot <- T1_beta_plot + theme_classic()+ guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75)
T1_beta_plot
Warning: Removed 5 rows containing non-finite values (stat_boxplot).
Warning: Removed 5 rows containing missing values (geom_point).
my_legend <- get_legend(T1_beta_plot)
Warning: Removed 5 rows containing non-finite values (stat_boxplot).
Warning: Removed 5 rows containing missing values (geom_point).
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_beta_legend.pdf")
Saving 7.29 x 4.51 in image
T1_beta_plot<-T1_beta_plot+ theme(legend.position = "none")
T1_beta_plot
Warning: Removed 5 rows containing non-finite values (stat_boxplot).
Warning: Removed 5 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T1_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 5 rows containing non-finite values (stat_boxplot).
Warning: Removed 5 rows containing missing values (geom_point).
T1_beta_df<- T1_beta$data
T1_betamod<-aov(formula = beta_div_value ~ group ,data = T1_beta_df)
summary(T1_betamod)
Df Sum Sq Mean Sq F value Pr(>F)
group 4 0.2967 0.07419 8.83 1.1e-06 ***
Residuals 250 2.1003 0.00840
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(x = T1_betamod, which = "group")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ group, data = T1_beta_df)
$group
diff lwr upr p adj
Dicamba-Atrazine-Mesotrione -0.027993939 -0.07861966 0.022631781 0.5508308
Glyphosate-Atrazine-Mesotrione -0.036412121 -0.08703784 0.014213599 0.2806138
Handweeded-Atrazine-Mesotrione -0.058666667 -0.11176337 -0.005569963 0.0220819
Non-Treated-Atrazine-Mesotrione 0.040006061 -0.01061966 0.090631781 0.1940292
Glyphosate-Dicamba -0.008418182 -0.05644596 0.039609594 0.9889807
Handweeded-Dicamba -0.030672727 -0.08129845 0.019952993 0.4576855
Non-Treated-Dicamba 0.068000000 0.01997222 0.116027775 0.0012035
Handweeded-Glyphosate -0.022254545 -0.07288027 0.028371175 0.7468226
Non-Treated-Glyphosate 0.076418182 0.02839041 0.124445957 0.0001747
Non-Treated-Handweeded 0.098672727 0.04804701 0.149298448 0.0000019
T2_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Herbicide")
T2_beta_plot <- T2_beta$plot
T2_beta_plot <- T2_beta_plot+ theme_classic() + theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75)
T2_beta_plot
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T2_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
T2_beta_df<- T2_beta$data
T2_betamod<-aov(formula = beta_div_value ~ group ,data = T2_beta_df)
summary(T2_betamod)
Df Sum Sq Mean Sq F value Pr(>F)
group 4 0.080 0.019995 4.367 0.00201 **
Residuals 231 1.058 0.004578
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(x = T2_betamod, which = "group")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ group, data = T2_beta_df)
$group
diff lwr upr p adj
Dicamba-Atrazine-Mesotrione -0.029188889 -0.070787897 0.0124101197 0.3048204
Glyphosate-Atrazine-Mesotrione -0.020375758 -0.057770485 0.0170189703 0.5647456
Handweeded-Atrazine-Mesotrione 0.022206061 -0.015188667 0.0596007884 0.4780643
Non-Treated-Atrazine-Mesotrione -0.015555556 -0.054775477 0.0236643659 0.8113359
Glyphosate-Dicamba 0.008813131 -0.031069709 0.0486959716 0.9738007
Handweeded-Dicamba 0.051394949 0.011512109 0.0912777898 0.0043255
Non-Treated-Dicamba 0.013633333 -0.027965675 0.0552323419 0.8962631
Handweeded-Glyphosate 0.042581818 0.007106064 0.0780575719 0.0097763
Non-Treated-Glyphosate 0.004820202 -0.032574526 0.0422149299 0.9966025
Non-Treated-Handweeded -0.037761616 -0.075156344 -0.0003668883 0.0465006
T3_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Herbicide")
T3_beta$plot #+ scale_color_manual(values = c("#F8766D", "#A3A500", "#00BF7D", "#00B0F6", "#E76BF3")) +
T3_beta_plot <- T3_beta$plot
T3_beta_plot <- T3_beta_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("")
T3_beta_plot <-T3_beta_plot + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T3_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
T3_beta_df<- T3_beta$data
T3_betamod<-aov(formula = beta_div_value ~ group ,data = T3_beta_df)
summary(T3_betamod)
Df Sum Sq Mean Sq F value Pr(>F)
group 4 0.0672 0.016792 1.689 0.153
Residuals 252 2.5054 0.009942
TukeyHSD(x = T3_betamod, which = "group")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ group, data = T3_beta_df)
$group
diff lwr upr p adj
Dicamba-Atrazine-Mesotrione 0.009452020 -0.04731207 0.06621612 0.9909320
Glyphosate-Atrazine-Mesotrione -0.037518182 -0.08753732 0.01250096 0.2404830
Handweeded-Atrazine-Mesotrione -0.014736364 -0.06475550 0.03528278 0.9275528
Non-Treated-Atrazine-Mesotrione -0.023409091 -0.07637301 0.02955483 0.7430333
Glyphosate-Dicamba -0.046970202 -0.10570358 0.01176317 0.1840580
Handweeded-Dicamba -0.024188384 -0.08292176 0.03454499 0.7896781
Non-Treated-Dicamba -0.032861111 -0.09412180 0.02839957 0.5804264
Handweeded-Glyphosate 0.022781818 -0.02946147 0.07502511 0.7524585
Non-Treated-Glyphosate 0.014109091 -0.04096017 0.06917835 0.9554583
Non-Treated-Handweeded -0.008672727 -0.06374199 0.04639653 0.9926689
library(ggpubr)
ggarrange(T1_beta_plot, T2_beta_plot, T3_beta_plot, ncol = 1)
Warning: Removed 5 rows containing non-finite values (stat_boxplot).
Warning: Removed 5 rows containing missing values (geom_point).
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_within_group_beta.pdf", width = 5, height = 10)
Examination of dissimliarity across time points by treatment and then again by all chemical treatments combined.
T1_beta_df$Time<-"T1"
T2_beta_df$Time<-"T2"
T3_beta_df$Time<-"T3"
beta_div_T1_T2_T3 <- rbind(T1_beta_df, T2_beta_df, T3_beta_df)
beta_anova<-function(data, Herbicide = "Herbicide"){
df_sub<- data %>% filter(group == Herbicide)
mod<-aov(beta_div_value ~ Time, data = df_sub)
print(summary(mod))
print(TukeyHSD(mod, "Time"))
boxplot(df_sub$beta_div_value ~ df_sub$Time)
}
beta_anova(beta_div_T1_T2_T3, Herbicide = "Non-Treated")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.234 0.11699 13.48 4.36e-06 ***
Residuals 142 1.232 0.00868
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 -0.09016162 -0.13451209 -0.04581115 0.0000111
T3-T1 -0.07307273 -0.11742320 -0.02872226 0.0004285
T3-T2 0.01708889 -0.02942628 0.06360405 0.6599624
beta_anova(beta_div_T1_T2_T3, Herbicide = "Handweeded")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.0558 0.027878 4.704 0.0104 *
Residuals 152 0.9009 0.005927
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 0.04627273 0.009645196 0.08290026 0.0090758
T3-T1 0.03427273 -0.002354804 0.07090026 0.0719037
T3-T2 -0.01200000 -0.046747927 0.02274793 0.6928413
beta_anova(beta_div_T1_T2_T3, Herbicide = "Dicamba")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.0733 0.03664 5.34 0.00596 **
Residuals 124 0.8507 0.00686
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 -0.03579495 -0.07791783 0.006327933 0.1126414
T3-T1 0.02778838 -0.01433450 0.069911266 0.2647095
T3-T2 0.06358333 0.01727131 0.109895354 0.0041226
beta_anova(beta_div_T1_T2_T3, Herbicide = "Atrazine-Mesotrione")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.0292 0.01458 1.3 0.276
Residuals 153 1.7163 0.01122
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 -0.034600000 -0.08744615 0.01824615 0.2708607
T3-T1 -0.009657576 -0.05811807 0.03880292 0.8847098
T3-T2 0.024942424 -0.02351807 0.07340292 0.4442906
beta_anova(beta_div_T1_T2_T3, Herbicide = "Glyphosate")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.0096 0.004779 0.804 0.449
Residuals 162 0.9632 0.005945
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 -0.01856364 -0.05334455 0.01621728 0.4184568
T3-T1 -0.01076364 -0.04554455 0.02401728 0.7448495
T3-T2 0.00780000 -0.02698091 0.04258091 0.8564946
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide
index <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Handweeded", "Non-Treated")
sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])
#+ scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T1_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Mode")
T1_beta_chemical_combined_plot <- T1_beta_chemical_combined$plot
T1_beta_chemical_combined_plot<- T1_beta_chemical_combined_plot + theme_classic() + guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T1_beta_chemical_combined_plot
Warning: Removed 353 rows containing non-finite values (stat_boxplot).
Warning: Removed 354 rows containing missing values (geom_point).
my_legend <- get_legend(T1_beta_chemical_combined_plot)
Warning: Removed 353 rows containing non-finite values (stat_boxplot).
Warning: Removed 353 rows containing missing values (geom_point).
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_beta_combined_legend.pdf")
Saving 7.29 x 4.51 in image
T1_beta_chemical_combined_plot<-T1_beta_chemical_combined_plot+ theme(legend.position = "none")
T1_beta_chemical_combined_plot
Warning: Removed 353 rows containing non-finite values (stat_boxplot).
Warning: Removed 353 rows containing missing values (geom_point).
T2_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Mode")
T2_beta_chemical_combined_plot <- T2_beta_chemical_combined$plot
T2_beta_chemical_combined_plot<- T2_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T2_beta_chemical_combined_plot
Warning: Removed 373 rows containing non-finite values (stat_boxplot).
Warning: Removed 374 rows containing missing values (geom_point).
T3_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Mode")
T3_beta_chemical_combined_plot <- T3_beta_chemical_combined$plot
T3_beta_chemical_combined_plot<- T3_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T3_beta_chemical_combined_plot
Warning: Removed 348 rows containing non-finite values (stat_boxplot).
Warning: Removed 349 rows containing missing values (geom_point).
ggarrange(T1_beta_chemical_combined_plot, T2_beta_chemical_combined_plot, T3_beta_chemical_combined_plot, ncol = 1)
Warning: Removed 353 rows containing non-finite values (stat_boxplot).
Warning: Removed 354 rows containing missing values (geom_point).
Warning: Removed 373 rows containing non-finite values (stat_boxplot).
Warning: Removed 375 rows containing missing values (geom_point).
Warning: Removed 348 rows containing non-finite values (stat_boxplot).
Warning: Removed 350 rows containing missing values (geom_point).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_within_group_beta_chemical_combined.pdf", width = 5, height = 10)
T1_beta_df_chemical_combined <- T1_beta_chemical_combined$data
T2_beta_df_chemical_combined<- T2_beta_chemical_combined$data
T3_beta_df_chemical_combined<- T3_beta_chemical_combined$data
T1_beta_df_chemical_combined$Time<-"T1"
T2_beta_df_chemical_combined$Time<-"T2"
T3_beta_df_chemical_combined$Time<-"T3"
m1<-aov(beta_div_value ~ group, data = T1_beta_df_chemical_combined)
summary(m1)
TukeyHSD(m1, "group")
boxplot(beta_div_value ~ group, data = T1_beta_df_chemical_combined)
m2<-aov(beta_div_value ~ group, data = T2_beta_df_chemical_combined)
summary(m2)
TukeyHSD(m2, "group")
boxplot(beta_div_value ~ group, data = T2_beta_df_chemical_combined)
m3<-aov(beta_div_value ~ group, data = T3_beta_df_chemical_combined)
summary(m3)
TukeyHSD(m3, "group")
boxplot(beta_div_value ~ group, data = T3_beta_df_chemical_combined)
beta_div__chemical_combined_T1_T2_T3 <- rbind(T1_beta_df_chemical_combined, T2_beta_df_chemical_combined, T3_beta_df_chemical_combined)
beta_anova(beta_div__chemical_combined_T1_T2_T3, Herbicide = "Chemical")
beta_anova(beta_div__chemical_combined_T1_T2_T3, Herbicide = "Hand")
beta_anova(beta_div__chemical_combined_T1_T2_T3, Herbicide = "Non-Treated")
treatment to control
plotDistances = function(p, m, s, d) {
# calc distances
wu = phyloseq::distance(p, m)
wu.m = melt(as.matrix(wu))
# remove self-comparisons
wu.m = wu.m %>%
filter(as.character(Var1) != as.character(Var2)) %>%
mutate_if(is.factor,as.character)
# get sample data (S4 error OK and expected)
sd = data.frame(sample_data(p)) %>%
select(s, d) %>%
mutate_if(is.factor,as.character)
sd$Herbicide <- factor(sd$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
# combined distances with sample data
colnames(sd) = c("Var1", "Type1")
wu.sd = left_join(wu.m, sd, by = "Var1")
colnames(sd) = c("Var2", "Type2")
wu.sd = left_join(wu.sd, sd, by = "Var2")
#remove this line to plot all comparisons.
wu.sd = wu.sd %>% filter(Type1 == "Hand" | Type1 == "Non-Treated")
# plot
ggplot(wu.sd, aes(x = Type2, y = value)) +
theme_bw() +
geom_point() +
geom_boxplot(aes(color = ifelse(Type1 == Type2, "red", "black"))) +
scale_color_identity() +
facet_wrap(~ Type1, scales = "free_x") +
theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
ggtitle(paste0("Distance Metric = ", m))
}
a<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T1"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
a <- a + ggtitle("Time 1 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T1_rare_allgroup_beta.pdf")
b<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T2"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
b <-b + ggtitle("Time 2 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T2_rare_allgroup_beta.pdf")
c<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T3"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
c<- c + ggtitle("Time 3 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T3_rare_allgroup_beta.pdf")
library(ggpubr)
ggarrange(a, b, c, ncol = 1)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_allgroup_beta.pdf", width = 7, height = 12)
Taxon abundance bar plot
#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
"#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
"#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
"#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
"#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
"#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
"#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
"#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
"#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
"#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
"#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
"#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
"#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
"#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
"#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
"#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
"#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
"#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
"#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
"#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
"#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
"#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
"#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
"#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",
"#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
"#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
"#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
"#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
"#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
"#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
"#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
"#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",
"#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
"#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")
phylumGlommed <- tax_glom(ps_rare, "Phylum")
#t1
phylumGlommed_herb_t1 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T1"), group = "Herbicide")
phylumGlommed_herb_t1 <- transform_sample_counts(phylumGlommed_herb_t1, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t1)$Herbicide <- factor(sample_data(phylumGlommed_herb_t1)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t1, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_Taxon_barplot_t1.pdf")
#t2
phylumGlommed_herb_t2 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T2"), group = "Herbicide")
phylumGlommed_herb_t2 <- transform_sample_counts(phylumGlommed_herb_t2, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t2)$Herbicide <- factor(sample_data(phylumGlommed_herb_t2)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t2, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/ITS_Taxon_barplot_t2.pdf")
#t3
phylumGlommed_herb_t3 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T3"), group = "Herbicide")
phylumGlommed_herb_t3 <- transform_sample_counts(phylumGlommed_herb_t3, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t3)$Herbicide <- factor(sample_data(phylumGlommed_herb_t3)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t3, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/ITS_Taxon_barplot_t3.pdf")
Combined herbicide and time bar plot
sample_data(ps_rare_sub)$herb_time<-paste(sample_data(ps_rare_sub)$Herbicide, sample_data(ps_rare_sub)$Time, sep = "_")
ps_rare_for_barplot <- prune_taxa(taxa_sums(ps_rare_sub) > 50, ps_rare_sub)
plot_bar(ps_rare_for_barplot, x= "herb_time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_BarPlot_Herbicide_Time.pdf", width = 20, height = 11)
Linear modeling of abundant taxa
Tax_glom_Subset <- function (physeq, y = "taxLevel", nreturns = "Number of returns"){
ps_1<- tax_glom(ps_rare_sub, taxrank = y )
myTaxa <- names(sort(taxa_sums(ps_1), decreasing = TRUE)[1:nreturns])
ps_1_sub <- prune_taxa(myTaxa, ps_1)
return(ps_1_sub)
}
ps_rare_family_top25<-Tax_glom_Subset(physeq = ps_rare, nreturns = 25, y = "Family")
ps_rare_order_top10<-Tax_glom_Subset(physeq = ps_rare, nreturns = 10, y = "Order")
#explore top 25 taxa with plot bar
plot_bar(ps_rare_family_top25, x= "herb_time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
plot_bar(ps_rare_order_top10, x= "herb_time", fill = "Order") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
#write function to wrangle data prior to anova
abund_aov_wrangle <- function (physeq, y = "Tax_Level"){
tax<-tax_table(physeq)[,y]
meta<-data.frame(sample_data(physeq))
counts<-data.frame(otu_table(physeq))
rownames(counts) <- tax[,1]
counts<-data.frame(t(counts))
counts$Time <- meta$Time
counts$Herbicide <- meta$Herbicide
counts$Herb_time <- meta$herb_time
return(counts)
}
test<-abund_aov_wrangle(ps_rare_family_top25, y = "Family")
mod_abund<-function(count_tab, IV = "Groups to be tested") {
for(j in 1:length(unique(count_tab[,"Herbicide"]))){
data <- count_tab %>% filter(Herbicide == unique(count_tab$Herbicide)[j])
#change this to the number of returns from the tax_glom_subset function
for (i in 1:25) {
mod <- aov(unlist(data[i]) ~ matrix(data[,IV]))
#sanity check
#print(c(j, i))
if(summary(mod)[[1]][["Pr(>F)"]][1] <= 0.05) {
print(summary(mod))
print(c(as.character(unique(count_tab[,"Herbicide"]))[j], names(data)[i]))
boxplot(unlist(data[i]) ~ unlist(data[IV]), main =paste(names(data[i]), as.character(unique(count_tab[,"Herbicide"]))[j]), xlab= "Time", ylab="Abundance")
}
}
}
}
mod_abund(test, IV = "Time")
#explore
ps_rare_certabosidiacae <- subset_taxa(ps_rare, Family == "f:Ceratobasidiaceae")
plot_bar(ps_rare_certabosidiacae, x= "herb_time") + scale_fill_manual(values = col_vector)